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This paper presents a novel framework for the modeling of biological networks. It makes use of 
recent tools analyzing the robust satisfaction of properties of (hybrid) dynamical systems. The main 
challenge of this approach as applied to biological systems is to get access to the relevant parameter 
sets despite gaps in the available knowledge. An initial estimate of useful parameters was sought by 
formalizing the known behavior of the biological network in the STL logic using the tool Breach. 
Then, once a set of parameter values consistent with known biological properties was found, we tried 
to locally expand it into the largest possible valid region. We applied this methodology in an effort to 
model and better understand the complex network regulating iron homeostasis in mammalian cells. 
This system plays an important role in many biological functions, including erythropoiesis, resistance 
against infections, and proliferation of cancer cells. 

1 Introduction 

In cellular biology, the value of the parameters are uncertain (when a measurement does exist) or not 
known at all (unmeasured or not computable from first principles). It is well known that values of kinetic 
constants obtained from in vitro measurements on purified enzymes can be significantly different from 
the in vivo values, due to interactions with other cell components, to regulatory modifications of enzymes, 
to sequestration, to anomalous diffusion or still other phenomena. In addition, the values available from 
the literature are heterogeneous (measurement performed in different species, or on different cell types, or 
in different conditions). Two sets of measurements performed on identical cell types placed in supposedly 
identical conditions can also be qualitatively different for multiple reasons (undetected heterogeneity 
of cell populations, different batches of antibodies used for detection, for instance). The fact that the 
available data and pieces of information are generally too scarce with respect to the size of dynamical cell 
models forces the modeler to include every bit of available data, at the cost of increased heterogeneity 
and thus increased uncertainty. But even with such an inclusive strategy the amount of data remains 
inadequate to identify a fully instantiated model. 

The complexities and heterogeneity of both the parameter values and the analytical expressions of the 
terms entering the differential equations describing a biological system implies such a system results from 
a family of dynamical systems, instead of a single one. The properties which formalize experimentally 
observed behaviors should thus be robust with respect to the uncertainties mentioned above. 
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These issues are considered in the present paper and the semi-formal method implemented to explore 
the parameter space is applied to an important question in biology: how iron homeostasis is maintained 
in mammalian cells? Iron is the most abundant transition metal required by most forms of life. The 
largest part needed by animals ends up in red blood cells, but all cell types do use the metal J2j. Iron 
is handled by specific molecules in animal physiology, so as to avoid the deleterious reactions the iron 
ions catalyze, conversion of oxygen derivatives into reactive species (the reactive oxygen species, ROS) 
prominent among them The mandatory requirement for iron and the negative impact the metal may 
have on cellular function call for a tight regulation of its use Q. Deregulation translates into cellular 
oxidative damage: this links iron homeostasis to the redox balance in animal, and other, cells. These 
connections underlie scores of (human) pathological situations in which these homeostatic systems (iron 
and redox homeostasis) are perturbed. 

Our approach rests on the following points: we deal explicitly with the parameter space, and we do 
not focus on a single instantiation of the model parameters. When the parameter values are uncertain 
due to the fact that they have been measured in different cell types or in different conditions, the data are 
formally represented by inequalities. In other words, we do not assign the measured value directly to the 
parameter but rather impose inequality constraints on that parameter. If there is no measurement at all, 
bounds are chosen to define a physiologically reasonable domain. Observed behavioral properties can be 
formalized using a temporal logic formalism. Then the behaviors of the parametrized system are explored 
by sampling predefined volumes of parameter space (including initial conditions and model parameters), 
and the instantiated models that satisfy the temporal logic formula are identified. As illustrated below 
this very basic scheme needs some elaboration to be applied to the modeling of a complex, yet simplified, 
biological phenomenon. 

The biological network we have built for iron homeostasis is presented in section [2| The system of 
differential equations is introduced in section [3] and its dynamical properties explored in the following 
section. 

2 Biological network 

We are considering regulation of iron handling in cell types which acquire iron by endocytosis of the 
transferrin-transferrin receptor complex (most of them) and express an iron exporter in the form of fer- 
roportin. Proper supply of iron to animal cells is required to avoid deleterious effects of too much or too 
little iron as such imbalance jeopardizes the proper function of cells and may lead to death. Regulation 
of iron supply occurs at two levels, a systemic one and a local one organized around Iron Regulatory 
Proteins. The latter are the main focus of this study since they directly act on the molecular nodes of the 
considered network. 

2.1 Species 

The considered network is shown in Figure [T] A. It contains five species: Iron Regulatory Protein (IRP, 
no difference is made between the two known proteins displaying this function), iron (Fe), ferritin (Ft), 
the transferrin receptor (TfRl) and ferroportin (FPNla). Moreover, transferrin (Tf) is a protein which 
carries iron in the blood stream throughout the body, and which can bind zero, one, or two iron atoms. It 
is the main provider of iron to most cells and, as such, it is an input of the network. The circle labeled by 
Fe_out is the iron exported out of the cell. The circle labeled by Fe_cons represents the iron consumed 
for the cell needs. 
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A) B) 

Figure 1: A) Schematic representation of the main biological processes involved in the cellular homeo- 
static control of iron. 

The dashed arrows represent translation of mRNA into proteins. The arrow ending in an unfilled white 
arrow represents amplification of translation, while the lines with a broken ending arrowhead represent 
inhibition of translation. The lines ending with a combined perpendicular stroke and arrow represent iron 
transport through membranes. The arrow leading to a pink circle indicates degradation. Finally, the two 
regular arrows represent the loading/unloading of iron into/from the ferritins and the consumption of iron 
for the cell needs. The green rounded rectangles represent proteins, the green parallelograms represent 
mRNA, and the blue circles represent atoms. The yellow concave hexagon represents the transferrin 
receptor. This diagram was drawn with the software CellDesigner liT4l . 

B) Interaction graph defined by the equations. The signs on the arrows indicate whether the interaction 
is positive or negative. 

The IRP are proteins which can bind to a specific hairpin-like structure on mRNA, called Iron Re- 
sponsive Element (IRE). Depending on the position of the IRE on the non-coding sequence of mRNA, 
two main outputs may result: 

• If the IRE is located at the beginning of the mRNA, in the 5' untranslated region (UTR), the 
binding of the IRP prevents the translation machinery from binding to the mRNA and translating 
it, so the binding of the IRP leads to a reduced quantity of the corresponding protein; 

• If the IRE is located at the end of the mRNA (3' -UTR), the binding of the IRP inhibits the degra- 
dation of this mRNA, leading to an increase in the quantity of translated protein. 

In the presence of enough iron, the regulatory activity of the IRP is decreased by different molecular 
mechanisms [24]. 

The iron species (Fe) represents the cellular iron which is available for cellular processes. The actual 
form of this iron in biological cells is still a topic of discussion. A virtual species, which we call "Fe", is 
introduced in order to model the usable iron in the cell. If this amount of intracellular iron becomes too 
low, the cell cannot maintain vital functions and dies. 
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Ferritin is a protein which stores iron [3]. The complex mechanisms of iron loading into and release 
from ferritin are not considered here. The mRNA of both ferritin subunits contain an IRE on their 5'-UTR 
regions, and the IRP inhibit their translation. 

Transferrin receptors are located on the surface of the cell CD. The fixation of transferrin on its 
receptor leads to the endocytosis of the complex. This process involves several steps and releases iron 
into the cell. Transferrin receptor mRNA contain five IRE in the 3'-UTR region, which is stabilized upon 
IRP binding. The density of transferrin receptor at the cell surface is thus correlated with the IRP activity. 

Finally, ferroportin (FPNla) is an iron exporter located at the surface of some cells ll20l . A form of 
FPNla contains an IRE on its 5'-UTR mRNA region, so the translation is inhibited by the presence of 
IRP. 

2.2 Behavior 

A qualitative description of the system has been obtained through a large body of biological experiments 
and we selectively summarize them here. 

If the amount of iron is sufficient (iron-replete situation), the cell is in a stationary state. Assuming 
that the iron provision stays constant, the concentration of each protein described in Figure [T] A does not 
vary. In this state, the IRP activity is expected to be relatively low. In the case of too much iron (which 
is not studied here), the saturation level of transferrin increases, the IRP activity reaches a minimum and 
inhibits translation of the iron importer (TfRl) to minimize cellular iron input. 

From the iron-replete stationary state, if the cells become iron-depleted, the inhibition of IRP activity 
by iron decreases. Thus, the increasing IRP activity leads to an increase in the transferrin receptor 
concentration and to a decrease in ferroportin concentration. This, in turn, leads to an augmentation in 
the amount of iron entering the cell and a reduction in the amount of iron leaving the cell. Moreover, the 
iron loaded in ferritin is released to convert stored iron into a form available for biosynthetic purposes. 
This release allows the cell to face a temporary reduction in the level of iron (typically because of a 
reduction in the iron input) for a certain amount of time. When ferritin is depleted, and if there is still 
no iron input mediated by the receptor TfRl, the cell runs out of iron and enters in a cellular death 
process. The IRP regulating system is effective as long as the iron level remains low. When the iron level 
regains an acceptable level again, the IRP activity decreases, and the amount of filled ferritin goes up, 
regenerating the stock of stored iron. 

3 ODE modeling 

In this section we describe the differential equations that we derived to specifically describe the evolution 
of the concentration of each species over time. Some processes in the iron homeostasis network are in 
fact a composition of many reactions, and the choice of an analytical expression to represent them is not 
straightforward. This is the case of the iron entry mediated by TfRl, which involves endocytosis, and of 
the loading of iron in ferritin, for example. 

3.1 Iron equation 

In the cell, when the level of iron becomes low, some iron-proteins release their iron content, by protein 
degradation for instance, which contributes to refill the iron pool; however, only iron release from ferritin 
is considered in our model without affecting the general validity of the reasoning. The equation defining 
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the iron concentration is: 

dFe dFt 

-^T = k Fe jnput ■ TfRl ■ Tf sat - n Ff-^-- kFe-export ■ Fe ■ FPNla - k Fe _ com Fe (1) 

The import of iron into the cell is proportional to the concentration of transferrin receptor and to 
the average amount of iron bound to transferrin (called transferrin saturation): Tf sat . The second term 
describes the storage or release of iron due to the synthesis or degradation of ferritin. This variation 
is equal to the number of iron atoms per ferritin molecule (aka np t ) times the variation in the ferritin 
concentration (we consider that (i) iron release occurs only when ferritin is degraded (ii) synthesized 
ferritins are very quicky filled with iron). The third term describes the export of iron out of the cell as 
proportional to the iron and FPNla concentration. The export parameter is kp e _ export . Finally, within 
the cell, iron is used for many purposes; this consumption of iron is represented by the last term in our 
equation (kp e cons ■ Fe). For the sake of simplicity, we consider this consumption term (which represents a 
whole set of cellular processes) to be proportional to the iron concentration. The export of iron by FPN1 
forms that are not regulated by IRP can be accounted for by this term. 



3.2 IRP equation 

The equation describing the IRP concentration is: 
dIRP 

= k IR p_ pm d - k Fe ^ mP ■ sig + (Fe, Fe ^ IRP ) ■ IRP - k IRPjdeg ■ IRP (2) 

We define a constant rate of production of IRP by kjRp _ pro d- The concentration of iron regulates 
the activity of IRP. The degradation of IRP is described by a constant basal term (kiRp_de g • IRP) and an 
iron-related term (kp e ^jRp ■ sig + (Fe,8p e ^jRp) IRP). Then, if the iron level is significantly below the 
threshold dFe^mp, the degradation rate is kiRp_d eg ■ IRP. If the iron concentration is significantly above 
this threshold, the degradation rate is (kpe^iRP + kmp_d eg ) • IRP, where kp e ^iRP is the parameter describing 
the inhibition of IRP by iron. 

The sigmoid function is defined such that: 

x n 

sig + (x,6) = — 

y ' x n + d n 

where x is a concentration variable, 6 is the threshold and n defines the steepness of the sigmoid. In our 
model, we use n = 30, which corresponds to a very steep sigmoid. 



3.3 Ferritin equation 

To model the temporal evolution of ferritin, we assume that each ferritin protein is rapidly filled with iron 
atoms, right after its synthesis is completed. In other words the amount of empty ferritin is negligible. 
The variable Ft thus represents the concentration of filled ferritin proteins within the cell. The equation 
describing the evolution of the concentration of ferritin is: 

dFt 

—j^- = kpuprod — kjRp^p t ■ sig (IRP, dmp^Ft) — kptjeg • Ft (3) 

We consider a basal production rate described by kpt^md- The second term describes the regulation by 
IRP. We use a sigmoidal regulation here, but we could consider other dynamics. If the activity of IRP 
is above the threshold Qirp^fu the production rate is lowered by kjRp^pt- It follows that kiRp^pt has to 
be lower or equal to kp tjpro d'. as IRP inhibits translation, this inhibition should never lead to a negative 
production rate. Finally, the third term {kpt_d eg ■ Ft) describes the spontaneous degradation of ferritin. 
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3.4 Ferroportin equation 

The equation describing the ferroportin concentration is very similar to that of ferritin. The equation is: 
AFPNla 

j- = kFPNJa_pmd ~ klRP^FPNla ' s 'g (IRP, QlRP^FPNla) ~ ^FPNlajleg ' FPNla (4) 

The first term describes the basal production of FPNla. The second term expresses the regulation 
of IRP on the translation of FPNla. With the same reasoning as previously, the regulation parameter 
kjRp^ppNia has to be lower or equal to the production parameter kppNia^md- The dynamics of IRP 
regulation on Ft and FPNla are considered similar because both FPNla and Ft mRNAs have a single 
IRE located in the 5'-UTR region without experimental evidence of major differences between the two. 



3.5 Transferrin receptor equation 

The concentration of the transferrin receptor evolves according to the following equation: 

— = krfRijprod. + kmp^TfRi ■ IRP — krfRijieg ■ TfRl (5) 

The production of TfRl includes a basal rate krfRijmd an d it is increased proportionally to the IRP 
concentration. This is taken into account in the term kjup^jfRi ■ IRP- This term is different from those 
describing the regulation of ferritin and ferroportin in an effort to consider the binding of several IRPs to 
a single TfR mRNA molecule (since it contains several IREs). The half-life of TfRl mRNA is considered 
proportional to the IRP concentration, although this is a rough approximation fl3l . The last term of the 
equation represents the spontaneous degradation. 



3.6 Interaction graph 

The interaction graph is a useful tool to visualize the pattern of interactions in a system of differential 
equations. It is defined from the Jacobian matrix of the system. Each non-zero element in this matrix 
indicates that a species influences another species. The nodes of the interaction graph are the dynamical 
variables of the system (concentration of biological species), and each arc represents the influence of a 
node on another, i.e. a non-zero element of the Jacobian matrix. The sign of this element determines 
whether the influence is positive or negative. Of course it is possible that the sign of the elements depends 
on the point of phase space at which the Jacobian is evaluated, but in our model the signs are constant. 
The interaction graph we obtain is shown in Figure [I] B. 

We define a negative (resp. positive) circuit (closed path) as a succession of arrows in the interaction 
graph such that the product of the signs labelling these arrows is negative (resp. positive). This graph 
contains three negative circuits of length greater than one (IRP— s-Fe— » IRP ; IRP— ^FPNla— »Fe— s-IRP and 
IRP— > TfRl— s-Fe— »IRP), and five negative loops corresponding to spontaneous degradation or consump- 
tion. It also contains one positive circuit ( IRP— >Ft— >Fe— >IRP ). The presence of these circuits, together 
with the presence of non-linear terms in the equations, generate a complex dynamical behavior. 



4 Methodology 

In this section we present our strategy to study the system. We first translate all available biological data 
into inequalities on parameters and variables. Then, we formally define the expected behavior for the 



two modes described in Section 2.2 the iron-replete mode and the iron-depleted mode. Finally we look 



for a subspace which robustly verifies this behavior. 
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4.1 Parameter search space 

The parameter space of the system is twenty-dimensional. To define an initial parameter search space, we 
combined published data from different sources as shown in table [T] Whenever possible, the available 
data obtained with erythroleukemic cells have been used, as these cells are relatively immature and 
proliferate, as considered for the initial steady state. 



In mouse macrophages, half-life of FPNla 
seems to be more than 20 hours lfT8l 


< k FPNla _ deg < 9.6 x 10~ 6 s- 1 


The production rate of TfRl was estimated at 
7x 10" 6 pg/cell/min as in [22] 


1.0 x 10~ 13 < k TfR1 _p rod < 2.0 x 10~ 13 mol/L/s 


The half-life of IRP is 12-15 hours in H1299 hu- 
man lung cancer cells fTT] 


1.28 x 10~ 5 < k mP je g < 1.6 x 10~ 5 s" 1 


In human erythroleukemia cell line K562, half- 
life of TfRl is 8 hours |23 


2.0 x 10~ 5 < k TfRljdeg < 3.0 x 10~ 5 s" 1 


The TfRl production rate enhanced by IRP is 
6x 1(T 5 pg/cell/min E2 


4.2 x 10~ 5 < kjRp^Tmi < 14.4 x 10~ 5 s" 1 


In human erythroleukemia cell line K562, 
iron input is up to 7.2±2.4 x 10 4 Fe 
atoms/cell/min ll23l 


2 x 10- 2 < k FeJnput < 3.9 x 10- 2 s- 1 


In mammalian cells, the number of iron atoms 
in ferritin up to 4500 (3l 


< n Ft < 4500 


Human normal saturation level of transferrin is 
withm 25% and 45% [21 ], we take 30% 


Tfsat = 0,3 (in iron-replete situation) 


In human erythroid cells, a fraction of the iron 
exporter (FPN1) is not IRP-regulated [6]. The 
amount of iron needed for biosynthesis ("con- 
sumed") is higher than the exported one. 


k F e.cons > k Fe export -FPNla (in iron-replete situ- 
ation) 


A higher limit for the iron pool in ery- 
throleukemic cells is set at 2 /xmol/L lTT2l 


Fe < 2 x 10~ 6 mol/L 


In rats liver, the IRP concentration in around 
0.11 pmol/mg of cytosolic protein under normal 
conditions fl4j 


3 x 10" 9 < IRP < 10.7 x 10~ 9 mol/L (in iron- 
replete situation) 


In human erythroleukemic (K562) cells, there 
are 140,000 transferrin receptors by cell ifTTTl 


1.0 x 10~ 8 < TfRl < 10.0 x 10" 8 mol/L (in 
iron-replete situation) 



Table 1: The left column shows the data collected in the literature; the right column describes the corre- 
sponding interval on parameters and variables, and the inequalities deduced. 



We then obtained intervals for eight parameters, including one in the iron-replete situation. Moreover, 
we derived an inequality linking k Fe cons , k Fe _ exp ort and FPNla. Finally, we found data on the value of 
three variables {Fe, IRP and TfRl) among which two are related to the iron-replete situation. 

In order to evaluate the search space for the twelve other parameters, we proceeded by analogy. For 
example, we had no information on the basal production rate of ferritin, ferroportin and IRP, but we 
knew an interval for the production rate of the transferrin receptor. We then assumed that the basal pro- 
duction rate of ferroportin was neither more than 1000 times higher than the highest possible transferrin 
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receptor basal production rate, nor 1000 times lower than the lowest. Concerning the thresholds param- 
eters, we simply assumed that their possible intervals were the same as for the corresponding protein 
concentrations. 

4.2 The iron-replete steady state 

The iron-replete situation is modeled by a steady state, so we do not consider evolutions of the variables 
over time. In this state, the value of the variable representing the iron concentration is expected to be sig- 
nificantly higher than the threshold dpe^mp- We can thus infer that the sigmoidal term sig + (Fe, dpe^iRp) 
is very close to one, and we replace it. The value of the variable representing the IRP activity is also 
expected to be largely below the thresholds Qmp^Ft and Qirp^fpniu- We then consider that the two 
sigmoidal functions sig + (IRP, Qmp^Ft) and sig + (IRP, Oirp^fpnio) are very close to zero. These simpli- 
fications contribute to reducing the complexity of the system. Note that this amounts to approximating 
the sigmoid by step functions, and consequently the system could be viewed as an hybrid system. 

We developed two complementary approaches for studying the system in the iron-replete steady 
state. The first consists in a mathematical analysis of the system. The second consists in performing 
deductions from the quantitative data. 

Looking for a steady state, we set the derivatives in the model to zero and solved the resulting 
algebraic system using the symbolic solver Maxima lfT9l . This way, we could prove the existence of one 
unique stationary point, given by Equations (|6]> to ( fT0| ). 



^ _ Tf sat • kppNla.deg ' kp e _input ' (krfRl.prod ' (klRPjdeg + kp e -^IRP) + kjRp^pfRl ' klRP_p m d) ^ 
(kp e 

-export ' kppNla.prod ~t~ kp e cons ■ kppNi ac ieg) • 

\klRP_deg + kFe^lRp) • k T f R i_ deg 

Ft = (7) 

kFtjieg 

FPNla = k FPNla.prod (g) 
kFPNXajdeg 

IRP = k IRP _ prod ^ 

klRP_deg + kpe^IRP 
jij-j^y _ {klRPJeg + kpe^IRp) • kjfR\_ pro d + kjRp^pfRi ■ kjRp pro d 
(klRP_deg + kpe^iRp) ■ kpfRijgg 



We then considered the stability of this stationary point by computing its Jacobian matrix since an 
unstable stationary point cannot obviously represent an observed cellular state. The obtained Jacobian 



matrix is shown in equation (111. 



f FPNla ■ kp e export kp e cons Ff sat ' kp e Jnput Fe ■ kp e export kptjdeg' n Ft N 

—kpfRl_deg klRP^TfRl 

—kFPNlaJeg 

-k Ft _deg 

V —kmpdeg — kFe^lRp) 

(11) 

As this matrix is upper triangular, its eigenvalues are equal to its diagonal elements. It follows that all 
eigenvalues are real and negative, which proves the stability of the stationary point. 
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Another way to study the steady state consists in using the specific data we have on the variables at 
the stationary state and propagate them on the parameter space. For this purpose, we use Realpaver |[T6ll . 
an interval solver. The solver uses the algebraic equations between parameters and variables, as well as 
defined intervals and constraints, as inputs for the calculation. The aim of Realpaver is to reduce the 
provided intervals and consequently reduce the parameter space. The results are the following: 



k F PNla.deg > l.Ox 10 13 S 



< F PNla_prod 



MRPjrod 



<9.6x 10~ n mol/L/s 



>3.84x 10~ 14 mol/L/s 



[l.Ox 
[l.Ox 



k Fe ^IRP 



<3.39x 10~ 4 s~ ! 
< 3.33 x 10~ 2 s- 1 



, l.Ox 10~ 10 ]) 
, l.Ox lO" 10 ]) 

1.0]) 



(initial interval: [0.0, 9.6 x 10" 6 ] ) 
(initial interval: 
(initial interval: 
(initial interval: [0.0, 1.0]) 
(initial interval: [l.Ox 10~ 9 , 
Realpaver returns intervals which are reduced by many orders of magnitude for four of the par ame- 
ters. These results prove that some subspaces of the parameter space do not contain valid parameter sets, 
and then, need not be considered when searching for an parameter set producing an adequate iron-replete 
steady state. This will reduce the search of the parameter space in the next step (analysis of the dynamical 
behavior). 

Besides, Realpaver makes no difference between variables and parameters considered in the algebraic 
system, they are all variables (unknowns). Consequently it also propagates the constraints from the 
parameters to the variables, and provides deductions on three variables in the iron-replete situation: 

• FPNla > 1.04 x 10~ 13 mol/L (initial interval: [1.0 x 10~ 13 , 1.0 x 10~ 5 ] ) 

• Fe > 3.5 x 10~ 8 mol/L (initial interval: [0.0, 2.0 x 10~ 6 ] ) 

• TfRl < 8.7 x 10~ 8 mol/L (initial interval: [1.0 x 10~ 8 , 1.0 x 10~ 7 ] ) 

These results, provided by Realpaver in less than one second on a Core 2 Duo 3 GHz with 8 Gb of 
RAM, can be considered as deductions made from the initial data. Indeed, the difficulty in finding a valid 
parameter set is that steady state values must be within their intervals. These two approaches facilitate 
the manual search of an initial parameter space. The analytical solution of the equations indicates, for 
each parameter, the variables impacted by a change in the parameter value. We use this information to 
define an order for checking variables, such that there are parameters allowing us to adjust a variable 
value but not changing the value of the already correct variables. Consequently, to find an initial valid 
parameter set, one has only to tune some parameters to fix the value of the first variable, then tune some 
other parameters to fix the value of the second variable, and so on. If the value of one variable cannot 
be set in its interval, one has to step back and change the value of parameters used to fix the value of 
the previous variable. The variable order we use is: IRP, TfRl (through the adjustment of kjRp^TjRi, 
k T fRUprod and k T jRi_deg), FPNla, Fe (through k Fe _ com , k Fe _ export and k FeJnput ) and finally Ft. Using this 
method, we found the following initial parameter set, denoted as 



kjj-Rl-prod 
kTfRLdeg 

n Ft 

klRP_prod 
k F e_export 
QlRP^ F PNla 
klRP^ F PNla 
klRP^TfRl 
k F e_cons 
klRP-^ F t 



1.7 X 10~ 13 

2.4 x 10~ 5 
400 

8.0 x 10~ 12 
300 

3.0 x 10~ 8 
5.0 x 10~ 13 
1.4 x 10~ 4 
3.0 x 10~ 4 
7 x 10 



mol/L/s 



mol/L/s 
L/mol/s 
mol/L 
mol/L/s 



-1 



k F PNla-deg 
k F e_input 
k F t_prod 
k F t_deg 
k Fe ^IRP 
Q F e^IRP 
klRP_deg 
k F PNla_prod 
QlRP-^ F t 



5.0 x 10 _t) 
3.0 x 10~ 2 
7 x 10~ n 
5.0 x 10~ 3 
l.Ox 10" 3 
1.5 x 10~ 7 

1.4 x 10~ 5 

2.5 x 10~ 12 
3.0 x 10~ 8 



,-i 



mol/L/s 



s 

s- 1 
mol/L 

s- 1 

mol/L/s 
mol/L 



-ii 



mol/L/s 
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It will be used as the initial parameter set for the exploration of the parameter space, as explained in 
the following section. 

4.3 Dynamics in iron-depleted situation 

We are now interested in the behavior of the network when the input of iron is cut. As a consequence, 
we need to compute the evolution of the system over time. For that, we use Breach [HHU, a tool based 
on Matlab/C++. This tool allows not only to simulate differential systems, but also to express temporal 
logic formula and to check whether a trajectory satisfies or not a given formula. 

We introduce briefly STL (Signal Temporal Logic). A more complete description can be found in (5|. 
Initially created for verifying the correctness of reactive programs, temporal logics provide languages to 
express properties over time-dependent phenomena for any kind of dynamical system. STL and its 
quantitative semantics [ 10] specializes into systems with continuous time and with real-valued states. 
An STL formula has the following syntax: 



where fx is an inequality constraint between expressions involving variables (written var [t] ), constants 
or derivatives of a variable over time (written ddt{var} [t] ). A formula is evaluated at each requested 
time (thereafter called "the current time"). A formula with a and statement is true if both q>\ and (fa are 
true at the current time. The ev_ [a , b] statement is true if q> holds at least once between the current time 
plus a and the current time plus b. If the interval is omitted, then a is and b is inf, that is, the formula is 
true if cp holds at least once, anytime after the current time. The statement alw_ [a,b] is true if q> holds 
between current time plus a and current time plus b. Again, if the interval is omitted then a is and b is 
inf meaning in that case that the formula must hold all the time. 

The experiment we want to simulate is the cut-off of iron entry for cells initially in the iron-replete 
situation. The total duration of the simulation is 48 hours, as the qualitative description of the iron- 
depleted situation tells us that significant changes can be seen in 24 hours. First we want the simulations 
to stabilize in a steady state corresponding to the iron-replete situation. Then, at time t=6 hours, we cut 
the entry of iron by setting Tf sat to zero. We expect the system to turn on the IRP regulatory system in 
order to maintain a correct iron level for at least ten hours. To identify parameter values for which the 
simulations satisfy all theses properties, we use temporal logic formula. 

We first want the system to be in a steady state corresponding to the iron-replete situation. To avoid an 
inappropriately long stabilization phase, we specify initial conditions close to the value of the stationary 
state. To ensure that the system reaches a steady state, we use the following temporal logic formula: 



<p = /l I (<pi) and (<P2) I ev_[a,b] (<p) | alw_[a,b] (<p) 



(12) 



(p S \ = abs(ddt{Fe> [t] / Fe [t] ) < 1.0e-4 

(p s2 = abs (ddt{Ft> [t] / Ft [t] ) < 1 . Oe-4 

<p 53 = abs(ddt{FPNla}[t] / FPNla[t]) < 1.0e-4 

<p S4 = abs(ddt{IRP}[t] / IRP[t]) < 1.0e-4 

<p 55 = abs(ddt{TfRl}[t] / Tf Rl [t] ) < 1.0e-4 



(13) 



(To make the reading easier, all formula dealing with the iron-replete steady state will begin by "<ps", 
formula related to parameters will begin by "<pp" and formula related to the behavior will begin by "<ps".) 
To check that the value of each variable at the steady state is in the intervals shown in table [TJ or deduced 
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in|4.2| we use the following formula: 



(pse = (3 . 5e-8 < Fe [t] ) 

(p sl = (3.0e-9 < IRP[t]) and (IRP[t] < 1.07e-8) 

<p s8 =(1.0e-8 < Tf Rl [t] ) and (Tf Rl [t] < 8.7e-8) (14) 

<p S9 = (1.04e-13 < FPNla[t]) and (FPNla[t] < 1.0e-5) 

(p sw = (1.0e-13 < Ft[t]) and (Ft [t] < 1.0e-5) 

As the maximal concentration of iron is 2 x 10~ 6 mol/L all the time, and not only at the steady state, we 
do not enforce this property here. 

Then, when the iron input is switched off, and the iron level becomes too low, the ferritins release 
iron. This release is due to a lower ferritin production rate, which provokes a decrease of ferritin con- 
centration. This can only be triggered by crossing the threshold dmp^Ft- This agrees with the previous 
statement that the variable describing IRP activity is lower than dmp^Ft under iron-replete conditions. 
We verify that property with the following formula: 

<p S ii=IRP[t] < theta_IRP_Ft (15) 

This increase of IRP activity can only be triggered by a decrease in iron concentration important enough 
to cross the threshold dFe^mp- This is in line with the already indicated inequality Fe > 0Fe^iRP under 
iron-replete conditions. The corresponding formula is: 

<Psi2 = Fe[t] > theta_Fe_IRP (16) 

Moreover, when the iron level becomes too low, the FPNla concentration is expected to become lower. 
This, also, can only be triggered by a crossing by IRP of the threshold Qirp^fpniq- That means that in 
the iron-replete steady state the FPNla concentration is above this threshold, which is expressed by this 
formula: 

<p 5 i3 = IRP[t] < theta_IRP_FPNla (17) 



Nevertheless , as explained in section |34| (resp. section [33| >, the regulation by IRP on FPNla (kjpp^ppNia) 
(resp. on Ft (kmp^pt)) cannot be higher that the basal production rate of FPNla (kFPNia-pmd) (resp. Ft 
(kFt_pmd))- The formula imposing that are: 

<p P i = k_IRP_FPNla <= k_FPNla_prod (18) 

and 

<p P2 = k_IRP_Ft <= k_Ft_prod (19) 

The activation of the regulatory system leads to the conservation of the iron level. We express that by 
enforcing the existence of a plateau of this level after shutting the iron import off. We want this plateau 
to exist for at least ten hours. To avoid incorrect simulations for which the concentration of this plateau 
is zero, we set it at, at least, one tenth of the concentration in the iron-replete situation. The following 
property expresses that: 

<Psi = ev_ [6*3600, inf] (alw_[0, 10*3600] 

( (<p S i) and (Fe[t] > . 01*Fe [4*3600] ) )) ( "°' 
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When the iron stored in ferritin is exhausted, the iron concentration eventually decreases to zero. 
Ferritin disappearance precedes iron depletion which is not rescued by replenishment from any internal 
store, in contrast to the transient phase immediately following removal of iron input from transferrin 
(Figure [2] A. Otherwise, it would mean that the cell has stored iron but is not using it, which contradicts 
the conditions set in |3.1| This implies that the parameter describing the strength of the regulation by the 
IRP (i.e. kjRp^.f t ) is slightly lower or equal to the production rate (i.e. kp t _ pm d)- We input this constraint 
with the following formula: 

<p P3 = k_Ft_prod*0.95 <= k_IRP_Ft (21) 

Next, in table [TJ two properties are expressed: the first enforces that the iron concentration is never 
higher than 2 x 10~ 6 mol/L, and the second enforces that the inequality kp econs > kp eexport ■ FPNla 
holds in iron-replete situation. The corresponding formula are: 

<p B2 = alw (Fe[t] < 2e-6) (22) 

and 

<p sl4 = k_Fe_cons > k_Fe_export*FPNla[t] (23) 

Finally, we have to define formula linking all the previous formula. For this purpose, we first define 
the formula tysall which aggregates all formula related to the steady state. As we want to be sure that, 
before the iron input is switched off, the system reaches and stays at the steady state, we enforce that 
the properties related to the steady state hold at least for one hour before the sixth hour. This formula 
expresses that: 

(pSall = ev_ [0,6*3600] (alw_ [0 , 3600] (<psi and ((p S2 and (<p S3 and ( ... and (psu)))) 

(24) 

We also define the formula (pBPaii enforcing that all formula related to parameters or behavior are satis- 
fied: 

<pBPall = <pPi and (<PP2 and ((p P3 and {(p B \ and <?>B2))) (25) 
And finally, the formula q> a u, which enforces that all properties are satisfied: 

<Pall = (<Psa//) and {(PBPaii) (26) 

It is noticeable that no information on TfRl in the iron-depleted situation nor on the parameters 
acting on it appears in this analysis. This is because, as we totally cut the entry of iron into the cell, 
variations in the concentration of TfRl do not change the amount of incoming iron. 



4.4 A robust valid region 

Having defined the formula which expresses the relevant properties, we search for a subspace where 
all parameter sets verify these properties, by first finding a valid instantiation of the parameters and 



then exploring a neighborhood of this instantiation. In section 4.2 we found an initial parameter set, 
satisfying constraints on the data for the steady state. It turns out that the simulation computed 
with this parameter set, shown in Figure [2] A, also verifies q> a u- As the temporal formula describes the 
expected qualitative behavior, this simulation is representative of the behavior of all trajectories in the 
region that we are looking for. We look around this point (in parameter space) to find a region exhibiting 
the expected behavior. To ensure that the formula q>P2 and <pp3 are satisfied, we enforce that kjpp^p t is 
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Figure 2: A) Characteristic simulation of the expected qualitative behavior. It shows that the system 
switches to a new mode when iron is lacking. The parameter set used is At time t=6 hours, Tf sat is 
set to zero. The iron concentration consequently decreases, leading to an increase of IRP activity. This 
increase leads to a decrease in Ft concentration, keeping the iron level constant. This also leads to a 
decrease in FPNla concentration and to an increase in TfRl concentration. During this phase when iron 
is kept constant, the concentration of iron is close to the threshold dFe^mp and IRP activity is close to 
the thresholds Qirp^fpniu and dmp^Ft- Once Ft concentration is zero, the iron concentration tends also 
to zero. 

B) Histogram showing the sensitivity of iron and ferritin concentration to variations of parameter value. 
These results are normalized: for example, a variation of 1% in kp e cons value leads to a variation of 0.2% 
in Fe value. To avoid incorrect sensitivity due do an initial transient phase of large amplitude, we do not 
consider the first three hours when computing the sensitivity. Moreover, as the behavior of the cell (i.e. 
death) changes when iron is totally exhausted, we consider the sensitivity only between t=3 hours and 
t=20 hours. 



computed equal to 0.97 x kp vro ^. Moreover, if the parameter kppNiajjmd is lower than kiRp^FPNia, we 
switch them to force that the formula q>p\ is satisfied. So, the dimension of the search space is slightly 
reduced. 

To guide the search, we compute the sensitivity of the system to parameters. We show in Figure|2]B 
the sensitivity of iron and ferritin concentration to a variation of parameters value. Sensitivity to ferro- 
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portin, transferrin receptor and IRP are also computed, but are not shown. We define the sensitivity to a 
parameter as the maximal sensitivity of the five variables of this parameter. This sensitivity is used in the 
following manner: if the system is very sensitive to a parameter, we will explore a close neighborhood 
of this parameter. On the opposite, if the system sensitivity to a parameter is low, we will explore a large 
neighborhood. Nevertheless, it has to be noted that the system can be very sensitive to a parameter, and 
that even with large variations on this parameter the qualitative property is still satisfied. This happens 
when the sensitivity which is detected affects an aspect of the trajectories which is not captured by the 
property of interest. 

To verify that all simulations corresponding to a parameter subspace exhibit the expected behavior, 
we performed intensive trajectory computations. We randomly picked 10000 parameter sets within this 
subspace and checked the truth value of the formula <p a y. By expanding as much as possible the parameter 
interval, we found an hyper-rectangle centered on &o such that all parameter sets in this region exhibit 
the expected behavior. Some parameters of could change up to 70% of their initial value, with an 
average of 34% for all, while conserving the behavior of the system. 

We observe that for the whole region of parameter space characterized in the iron-depleted situation, 
the FPNla concentration remains at a rather high level, implying the considered cell keeps its ability to 
export iron despite shortage of the metal. 

5 Conclusion 

Two models of the core iron homeostasis system have been published previously. Omholt et al ll22l aim 
at building a "unifying meta-model", that is a model which is not tied to a particular cell type. They 
include explicitly the two different Iron Regulatory Proteins, IRP1 and IRP2. Both proteins have the 
same mRNA targets, so they seem to be redundant, but they respond to the iron level in the cell through 
different molecular mechanisms. The model of Omholt et al l22l does not explicitly include ferritin nor 
ferroportin. In addition, all interactions in the network are switch-like (sigmoidal kinetics). A numerical 
value was assigned to all parameters except the sigmoid steepness, and the stability of the unique steady 
state was monitored as the sigmoid steepness is varied. The more recent paper of Chifman et al Q 
is closer to our approach. This model considers the same five biological species, but the analytical 
expression of their differential system is different from ours. The major differences are for the IRP, 
ferritin and ferroportin equations (note that the graph shown in figure 1 of Q does not follow the formal 
definition given here for an interaction graph). The mathematical analysis of the parametrized model 
of [5 ] is general in the sense that it does not rely on parameter instantiation. They showed that the model 
produces a single steady state, the stability of which was studied by performing simulations with sampled 
parameter values. 

Here we have chosen an approach intermediate between Omholt et al ll22ll and Chifman et al Q : a 
parametrized model has been built without imposing numerical values to the parameters. The available 
experimental data allowed us to apply inequality constraints between the parameters, together with be- 
havioral constraints expressed in a temporal logic formalism. We could thus study not only steady states, 
but also dynamical responses to perturbations, as, for instance, the sharp decrease of iron input, starting 
from an iron-replete cell state. The region of parameter space for which the model exhibits the observed 
behavior has been characterized. It appeared that, in the iron-depleted situation, the FPNla concentration 
remains at a rather high level, implying the considered cell kept its ability to export iron despite shortage 
of the metal. This is a significant outcome of modeling which justifies that another regulatory mechanism 
than that of the IRP may be needed under these conditions. This additional regulation may be provided 
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by FPNla degradation by hepcidin. 

We plan to extend the present model. For example the influence of iron on the inhibition of IRP 
activity involves the iron cluster biosynthesis pathway which will be introduced in the model. Also, the 
regulatory mechanisms of IRP could be described more precisely by integrating mRNA concentrations 
in the model. These extensions will be performed in a controlled and stepwise manner. The qualitative 
behavior will be compared with that of the core model presented here. This will give an insight into the 
robustness of the dynamical property of interest. 
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